% Matlab Code for replication of "Optimal Debt Maturity and Firm Investment"
% by J. Jungherr and I. Schott, September 2020
function [funval, output] = firm_optim_final(val,firmstate,PARMS,doall)
% Inputs:  Guess (K, Leverage)
% Take as given: the firm's state today 

debt        = firmstate.DEBT;  
zz          = firmstate.zz;  
W           = PARMS.W;  
SDF         = PARMS.SDF;  
Z           = PARMS.Z;  
TAU         = PARMS.TAU;
DELTA       = PARMS.DELTA;
C           = PARMS.C;
PSI         = PARMS.PSI;
ZETA        = PARMS.ZETA;
XI          = PARMS.XI; 
sigmaeps    = PARMS.sigmaeps;  
ETA         = PARMS.ETA;  
F           = PARMS.F;

if nargin <= 3
    doall           = 0;
end

GAMMA       = 1;  
kguess      = val(1);  
btildeguess = val(1)*val(2);  

alpha       = PSI*ZETA;
beta        = (1-PSI)*ZETA;
Lstar       = (beta*Z(zz)*kguess^alpha./W).^(1/(1-beta));

% Output tomorrow
Y           = Z(zz).*(kguess^PSI*Lstar.^(1-PSI)).^ZETA;

% Solve for EPSBAR 
epsguess    = ((btildeguess-kguess)/(1-TAU) - (Y - F - W.*Lstar - DELTA*kguess - C*btildeguess))/kguess;

% Define integral over epsilon
PHI_u       = 1/2 *(1+erf((epsguess-0)/(sigmaeps*sqrt(2)))); 

% Recovery value in case of default 
cond_mean_u =  - sigmaeps * ( 1/(sqrt(2*pi))*exp(-1/2*(epsguess/sigmaeps).^2)) ./ ...
    (1/2 *(1+erf(((epsguess)/(sigmaeps * sqrt(2))) )));
KTILDE      =  PHI_u.*kguess + (1-TAU).*(Y - F - W.*Lstar - DELTA*kguess).*PHI_u + ...
    (1-TAU)*kguess*PHI_u.*cond_mean_u ;
scrapvalue  = (1-XI)*KTILDE;

% break-even condition for issuance of short-term debt 
issuance_ltd =  (SDF.*((1-PHI_u)*(btildeguess-debt).*(GAMMA + C + 0) + ...
    (btildeguess-debt)*scrapvalue/(0+btildeguess)));

% Costs of issuing debt
costs       = ETA * (0 + max(0,btildeguess-debt))^2; 

cond_mean   = (1/(sqrt(2*pi))*exp(-1/2*((epsguess-0)/sigmaeps).^2))./(1 - ...
    (1/2 *(1+erf(((epsguess-0)/sigmaeps)/sqrt(2)))) )*sigmaeps;
QPRIME_val  = kguess*(1-TAU)*(1-PHI_u).*(cond_mean - epsguess);

% Compute short-term bond price 
P_ltd       = SDF.*((1-PHI_u).*(GAMMA+C + 0) + scrapvalue/(0+btildeguess));
if btildeguess == 0
    P_ltd = PARMS.p_fl;
end

% Compute etilde 
eprime      = kguess - issuance_ltd + costs;

funval      = - (-eprime  +                 (SDF.*QPRIME_val));

if eprime < 0
    funval =  99999; 
end

%%

if doall == 1
    % Define remaining outputs:
    output.P_ltd        = P_ltd;
    output.epsbar       = epsguess;
    output.ktilde       = KTILDE;
    output.eprime       = eprime;
    output.Lstar        = Lstar;
    output.V            = -funval;
    output.Kstar        = kguess;
    output.Levstar      = val(2);
    output.Lstar        = Lstar;
end

